This invention is in the field of seismic prospecting, and is more specifically directed to the analysis of seismic survey signals obtained in such exploration.
The use of seismic survey techniques is commonplace in the prospecting for subsurface reservoirs of hydrocarbons such as oil and gas. As is fundamental in the art, seismic surveys are typically obtained by imparting acoustic energy into the earth, either at a land surface or into the water in a marine survey, and detecting the energy at other surface or marine locations away from the source. The detectors sense the arrival of the acoustic energy after reflection from subsurface strata and interfaces, such that the time delay between the imparting of energy by the source and its receipt by the detectors will be indicative of the depth below the earth surface at which the reflecting interface is located. Conversion of the time-domain reflection signals into depth will thus provide a survey of the subsurface geology.
FIG. 1 illustrates a conventional land-based seismic survey in a portion of the earth having multiple subsurface strata 2, 4, 6, 8. In the example of FIG. 1, a reflective geological interface is present in the earth between subsurface strata 6,8 along which a particular depth point DP of interest for this description is located. Of course, many depth points, both along the interface between strata 6,8 and also from other interfaces in the region, will be analyzed in a typical survey.
In the typical land-based seismic survey, an array of receivers are deployed at the surface of the earth, and acoustic energy is imparted into the earth at varying surface source locations; in the typical marine seismic survey, the seismic source is towed by a vessel and activated at selected locations, and the receivers are either deployed in the water at fixed locations or are also towed by the source vessel. In either survey, the imparted energy is reflected from interfaces between subsurface strata, and this reflected energy is detected by the receivers in the array. In the simplified land-based example of FIG. 1, acoustic energy is reflected from depth point DP along multiple source-receiver paths that have varying angles of reflection at depth point DP. As such, the reflected energy is detected at receiver locations that have varying offset distances from the vertical, "zero-offset", path above depth point DP. The so-called "zero-offset" location is the surface location which is directly above depth point DP (i.e., the zero-offset location is the location at which reflected energy from depth point DP would reflect back to the source location). Typically, a zero-offset signal is not actually obtained in a survey, but is instead extrapolated from the other traces of varying offset.
Two source-receiver paths P.sub.n, P.sub.f are illustrated in the simplified example of FIG. 1. Near-offset path P.sub.n corresponds to energy imparted at source location S.sub.n and detected by receiver R.sub.n after being reflected from depth point DP, while far-offset path P.sub.f corresponds to energy imparted at source location S.sub.f and detected by receiver R.sub.f, also after reflection from depth point DP. The offset distances of paths P.sub.n, P.sub.f are illustrated in FIG. 1 as offsets x.sub.n, x.sub.f, respectively. Because the overall length of far-offset path P.sub.f is significantly greater than that of near-offset path P.sub.n, the travel time of acoustic energy along far-offset path P.sub.f is much greater than that along near-offset path P.sub.n, despite the paths having a common depth point DP.
A well-known technique in the art for eliminating random noise from the survey signals is referred to as common midpoint (CMP), or common depth point (CDP), gather and stack. In this technique, the signals from each source-receiver path having the same depth point are selected from the survey data, forming a group of traces known as a CMP (or, equivalently, CDP) gather. Elimination of random noise is then performed by summing the traces in the gather into a "stack". This summation will tend to reinforce the signal portion of the traces corresponding to the reflection event, which should be common among all traces having the same reflection midpoint, while the random noise will tend to cancel in the summation.
As noted above for the example of FIG. 1, in the CMP gather for depth point DP.sub.i, the seismic traces detected by receivers R.sub.n, R.sub.f will have different reflection times for energy reflected from the same reflection point DP, because their paths P.sub.n, P.sub.f are of different length. Specifically, the time at which energy reflected from depth point DP will appear in the trace for receiver R.sub.f will be later than the time at which energy reflected from depth point DP will appear in the trace for receiver R.sub.n. This time-delay effect of varying offset must be taken into account in the CMP stack process in order for the signal portions of traces of varying offset to properly align and provide an accurate indication of the depth of the reflector. This is typically handled by way of "normal moveout correction", or NMO (also NMOC), in which the traces corresponding to source-receiver paths of various lengths are time-shifted relative to one another so that their detected reflection events are aligned in time. The amount of the time shift for a given trace will, of course, depend upon its offset distance and upon the velocity with which the acoustic energy travels in the strata along the shot-receiver path. Normal moveout correction therefore requires the estimation or determination of a velocity, commonly referred to as the "stacking" velocity, for deriving the necessary time shift as a function of offset.
By way of further background, another well-known survey approach is commonly referred to as vertical seismic prospecting (VSP). In a VSP survey, energy is imparted into the earth from subsurface locations within a wellbore, and detected at surface detectors. The time delay between the imparting of energy by the wellbore source and its receipt by the surface detectors, in the VSP survey, will be indicative of attributes, such as acoustic velocities, of the subsurface strata, and may also be indicative of discontinuities and interfaces among the various strata. Of course, similarly as in the case of the surface-source seismic survey, the offset distance between the wellbore and the surface detectors will directly affect the time delay at which acoustic energy is detected at the surface. As such, normal moveout corrections are also necessary to process seismic signals obtained in VSP surveys. For purposes of clarity, the following discussion will refer to the conventional surface-based source seismic surveys, it being understood that similar concerns are applicable for moveout correction in VSP surveys.
In the ideal geological situation of a homogenous and isotropic medium, the wave surface of the acoustic energy imparted into the earth is spherical. Assuming reflection of the spherical wave surface from a horizontal reflecting surface in this ideal case, the relationship between arrival time t of the energy at a given surface detector R, and the offset distance x of detector R corresponds to a centered hyperbola; the relationship between offset distance x and travel time t is thus a straight line in the t.sup.2 versus x.sup.2 domain. The NMO time correction .DELTA.t.sub.x for a given seismic reflection event in a seismic trace may be derived according to the well-known Dix equation: ##EQU1## where t.sub.0 is the zero offset reflection time of a given reflection, where x is the offset distance of the trace to be corrected, and where V.sub.v is the acoustic velocity in the homogeneous isotropic medium. Solving this relationship for .DELTA.t.sub.x will produce the NMO time shift for the reflection event in the trace being corrected. In practice, the stacking velocity V.sub.v is typically determined as a root-mean-square velocity, or as a best-fit average velocity by way of a semblance analysis or other technique, to account for the multiple strata (e.g., strata 2, 4, 6 in the example of FIG. 1) through which the acoustic energy travels.
As noted above, the Dix equation hyperbola assumes the ideal condition of isotropic velocity in a homogenous medium. In practice, however, the subsurface of the earth often exhibits variations in velocity that depend upon the angle of propagation of the acoustic energy from the vertical, zero-offset, axis. In particular, the horizontal ray velocity (V.sub.h) generally differs from the vertical ray velocity (V.sub.v) along the zero-offset path. These effect of these variations, which is commonly referred to in the art as transverse isotropy, is due to the heterogeneity of velocities among the subsurface strata in the survey, and also to the anisotropic response of the earth to acoustic energy. Physical factors that cause anisotropic response include preferential orientation of rock properties, such as in sedimentary rocks with aligned crystals and sand grains, and also rocks with directional microcracks induced by regional stress fields. A simple modification of the Dix hyperbolic equation accounts for differences between the vertical and horizontal ray velocities V.sub.v, V.sub.h, respectively, as follows: ##EQU2## This expression of the Dix equation corresponds to an elliptical wave surface of acoustic energy traveling through the earth, with different velocity components in the vertical and horizontal directions.
The propagation of long-wavelength acoustic energy has been observed to exhibit apparent transverse isotropy in stratified regions of the earth, even in cases where each layer in the stratified region is itself individually isotropic. This apparent transverse isotropy is due to the dominance, at long offsets, of high velocity layers in the travel time of acoustic energy through multiple layers in the region, which causes the apparent horizontal velocity to be higher than the zero-offset, vertical, velocity. Sedimentary rocks that are formed through sand-shale stratification and other cyclic sedimentation processes typically exhibit transverse isotropy. In transversely isotropic regions, the accuracy of the Dix equation is practically limited to relatively short offset distances, such as offsets ranging up to on the order of the depth of the reflecting interface. For example, referring to FIG. 1, each of paths P.sub.n, P.sub.f pass through multiple strata 2, 4, 6. While the acoustic velocities are uniform within each of strata 2, 4, 6, the relative distances traversed by each path P.sub.n, P.sub.f vary with increasing offset distance. Therefore, even assuming that strata 2, 4, 6 are each individually transversely isotropic (i.e., having the same acoustic velocity in all horizontal directions therewithin), the stacking velocity will typically not be constant over the range of offsets obtained in the survey. Accordingly, the Dix equation hyperbolic assumption remains valid only over near-offset traces, as the time-distance relationship is no longer hyperbolic at large offsets in vertically anisotropic regions, such as those including transversely isotropic strata as shown in FIG. 1.
An example of the skewing of the horizontal velocity in transversely isotropic regions is illustrated in the plot of FIG. 2. FIG. 2 is a plot of vertical velocity V.sub.v versus horizontal velocity V.sub.h, illustrating the elliptical wavefront of acoustic energy, according to the Dix equation modified to incorporate different horizontal and vertical velocities V.sub.v, V.sub.h, as described above. In FIG. 2, ellipse E.sub.d is a plot of the wavefront of acoustic energy in accordance with the Dix equation, using the hyperbolic assumption described hereinabove. In this example, a horizontal velocity V.sub.he is estimated from the Dix equation. Skewed ellipse E.sub.s illustrates the actual horizontal velocity V.sub.ha that occurs in transversely isotropic media. As illustrated in FIG. 2, the actual horizontal velocity V.sub.ha is slightly larger than that estimated from the Dix equation according to the hyperbolic assumption. This discrepancy from the Dix equation will be most evident in far-offset traces, of course, considering the longer path distances presented therein. Accordingly, use of the Dix equation in conventional normal moveout processing is therefore not adequate to derive accurate normal moveout correction values for far-offset data in transversely isotropic regions, due to the discrepancy between the actual horizontal velocity and that which would be applied through solution of the Dix equation using conventional methods.
Conventional approaches for determining normal moveout corrections for far-offset seismic traces in transversely isotropic media are known. One approach utilizes a truncated Taylor series expansion to model the t.sup.2 -x.sup.2 relationship in the survey region. As described in Hake et al., "Three-term Taylor Series for t.sup.2 -x.sup.2 Curves of P- and S-Waves over Layered Transversely Isotropic Ground", Geophysical Prospecting 32 (1984), pp. 828-850, a three-term Taylor series expansion of t.sup.2 as a function of x.sup.2 may be used to calculate the elastic parameters of the transversely isotropic media. Quartic Taylor series expansion, using exact coefficients, has also been applied to the transversely isotropic layered media problem, as described in Tsvankin, et al., "Nonhyperbolic reflection moveout in anisotropic media", Geophysics, Vol. 59, No. 8 (August 1994), pp. 1290-1304.
The truncated Taylor series expansion has been successful in deriving nonhyperbolic moveout correction for long-offset reflection seismic data, particularly in improving the results of CMP stacking to improve the image of reflection signals. However, it is also desirable to derive estimates of the anisotropy parameters from the moveout analysis, to permit inversion of subsurface parameters that are indicative of the nature of the various strata in the survey region. It has been found, in connection with the present invention, that truncated Taylor series analysis is not well-suited for subsurface parameter inversion, because of the sensitivity of the Taylor series coefficients to measurement errors. The heightened error sensitivity is primarily due to the higher order (third order and higher) terms used in the Taylor series analysis.
By way of further background, another approach toward estimation of stiffness coefficients for multiple transversely isotropic layers in a seismic survey region is described in Byun et al., "Transversely isotropic velocity analysis for lithology discrimination", Geoexploration 28 (Elsevier Science Publishers B. V., 1991), pp. 269-292, incorporated herein by this reference. According to this approach, a skewed hyperbolic moveout formula is applied to pressure waves to derive three measurement parameters which, when combined with two shear wave measurement parameters, enable inversion for the five stiffness coefficients by way of a model-iterative optimization scheme.